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We present a new way of modeling turbulent thermonuclear deflagration fronts 
in Chandrasekhar-mass white dwarfs, consisting of carbon and oxygen, undergoing 
a type la supernova explosion. Our approach is a front capturing/tracking hybrid 
scheme, based on a level set method, which treats the front as a mathematical dis- 
continuity and allows for full coupling between the front geometry and the flow 
field. First results of the method applied to the problem of type la supernovae are 
discussed. It will be shown that even in 2-D and even with a physically motivated 
sub-grid model numerically "converged" results are difficult to obtain. 
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1. INTRODUCTION 

Type la supernovae, i.e. stellar explosions which do not have hydrogen in their spectra, 
but intermediate-mass elements, such as silicon, calcium, cobalt, and iron, have recently 
received considerable attention because it appears that they can be used as "standard can- 
dles" to measure cosmic distances out to billions of light years away from us. Moreover, 
observations of type la supernovae seem to indicate that we are living in a universe that 
started to accelerate its expansion when it was about half its present age. These conclu- 
sions rest primarily on phenomenological models which, however, lack proper theoretical 
understanding, mainly because the explosion process, initiated by thermonuclear fusion of 
carbon and oxygen into heavier elements are difficult to simulate even on supercomputers, 
for reasons we shall discuss in this article. 
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The most popular progenitor model for the average type la supernovae is a massive white 
dwarf, consisting of carbon and oxygen, which approaches the Chandrasekhar mass, Mch 
~ 1 .39 Mq, by a yet unknown mechanism, presumably accretion from a companion star, 
and is disrupted by a thermonuclear explosion (see, e.g., [52] for a review). Arguments in 
favour of this hypothesis include the ability of these models to fit the observed spectra and 
light curves. 

However, not only is the evolution of massive white dwarfs to explosion very uncertain, 
leaving room for some diversity in the allowed set of initial conditions (such as the tem- 
perature profile at ignition), but also the physics of thermonuclear burning in degenerate 
matter is complex and not well understood. The generally accepted scenario is that explo- 
sive carbon burning is ignited either at the center of the star or off-center in a couple of 
ignition spots, depending on the details of the previous evolution. After ignition, the flame 
is thought to propagate through the star as a sub-sonic deflagration wave which may or 
may not change into a detonation at low densities (around 10 7 g/cm 3 ). Numerical models 
with parameterized velocity of the burning front have been very successful, the prototype 
being the W7 model of [38]. However, these models do not solve the problem because 
attempts to determine the effective flame velocity from direct numerical simulations failed 
and gave velocities far too low for successful explosions [27, 19, 1], This has led to some 
speculations about ways to change the deflagration into a supersonic detonation [15, 16]. 

Numerical simulations of any kind of turbulent combustion have always been a chal- 
lenge, mainly because of the large range of length scales involved. In type la supernovae, 
in particular, the length scales of relevant physical processes range from 10~ 3 cm for the 
Kolmogorov-scale to several 10 7 cm for typical convective motions. In the currently fa- 
vored scenario the explosion starts as a deflagration near the center of the star. Rayleigh- 
Taylor unstable blobs of hot burnt material are thought to rise and to lead to shear-induced 
turbulence at their interface with the unburnt gas. This turbulence increases the effective 
surface area of the flamelets and, thereby, the rate of fuel consumption; the hope is that fi- 
nally a fast deflagration might result, in agreement with phenomenological models of type 
la explosions [38]. 

Despite considerable progress in the field of modeling turbulent combustion for astro- 
physical flows (see, e.g., [34]), the correct numerical representation of the thermonuclear 
deflagration front is still a weakness of the simulations. Methods used up to now are based 
on for the reactive-diffusive flame model [17], which artificially stretches the burning re- 
gion over several grid zones to ensure an isotropic flame propagation speed. However, the 
soft transition from fuel to ashes stabilizes the front against hydrodynamical instabilities 
on small length scales, which in turn results in an underestimation of the flame surface area 
and - consequently - of the total energy generation rate. Moreover, because nuclear fusion 
rates depend on temperature nearly exponentially, one cannot use the zone-averaged values 
of the temperature obtained this way to calculate the reaction kinetics. 

The front tracking method described here cures some of these weaknesses. It is based 
on the so-called level set technique which was originally introduced by Osher and Sethian 
[39]. They used the zero level set of a n-dimensional scalar function to represent (n — 1)- 
dimensional front geometries. Equations for the time evolution of such a level set which is 
passively advected by a flow field are given in [47] . The method has been extended to allow 
the tracking of fronts propagating normal to themselves, e.g. deflagrations and detonations 
[46]. In contrast to the artificial broadening of the flame in the reaction-diffusion-approach, 
this algorithm is able to treat the front as an exact hydrodynamical discontinuity. 
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In the following sections we present the main ideas and governing equations of this ap- 
proach. We describe only the most simple implementation of the flame model and leave ex- 
tensions including the reconstruction of the thermo-dynamical properties of mixed cells to 
subsequent papers. The main emphasis will be to demonstrate that the method can indeed 
be applied to the supernova problem. In addition, we will show that even if one attempts to 
model the physics of thermonuclear burning on unresolved scales well by physically mo- 
tivated LES one still has to perform calculations with very high spatial resolution in order 
to be at least near to a converged solution. 

The numerical results presented here were obtained by performing 2D simulations, but 
there is no problem extending them to 3D. Such simulations are under way. 

2. GENERAL PROPERTIES OF EXPLODING WHITE DWARFS 

To a very good approximation, the exploding white dwarf material can be described as 
a fully ionized plasma with varying degrees of electron degeneracy, satisfying the fluid ap- 
proximation. The governing equations are the hydrodynamical equations for mass, species, 
momentum, and energy transport including gravitational acceleration, viscosity, heat and 
mass diffusion [24], and nuclear energy generation [2]. They must be supplemented by an 
equation of state for an ideal gas of nuclei, an arbitrarily relativistic and degenerate electron 
gas, radiation, and electron-positron pair production and annihilation [8] . The gravitational 
potential is calculated with the help of the Poisson equation. In numerical simulations that 
fully resolve the relevant length scales for dissipation, diffusion, and nuclear burning it is 
possible to obtain the energy generation rate from a nuclear reaction network (for a recent 
overview, see [49]) and the diffusion coefficients from an evaluation of the kinetic trans- 
port mechanisms [32]. If, on the other hand, these scales are unresolved - as is usually the 
case in simulations on scales of the stellar radius - subgrid-scale models are required to 
compute (or parameterize) the effective large-scale transport coefficients and burning rates, 
which are more or less unrelated to the respective microphysical quantities [19, 34]. 

Initial conditions can be obtained from hydrostatic spherically symmetric models of the 
accreting white dwarf or - for Chandrasekhar mass progenitors - from the Chandrasekhar 
equation for a fully degenerate, zero temperature white dwarf [21]. Given the initial con- 
ditions and symmetries specifying the boundary conditions, the dynamics of the explosion 
can in principle be determined by numerically integrating the equations of motion. [31] 
gives a detailed account of some current numerical techniques used for modeling super- 
novae. 

2.1. Nuclear burning in degenerate C+O matter 

Owing to the strong temperature dependence of the nuclear reaction rates, S <~ T 12 
at T w 10 10 K [13], nuclear burning during the explosion is confined to microscopically 
thin layers that propagate either conductively as subsonic deflagrations ("flames") or by 
shock compression as supersonic detonations [7, 24]. Both modes are hydrodynamically 
unstable to spatial perturbations as can be shown by linear perturbation analysis. In the 
nonlinear regime, the burning fronts are either stabilized by forming a cellular structure or 
become fully turbulent - either way, the total burning rate increases as a result of flame 
surface growth [26, 50, 55]. Neither flames nor detonations can be resolved in explosion 
simulations on stellar scales and therefore have to be represented by numerical models. 

When the fuel exceeds a critical temperature T c where burning proceeds nearly instan- 
taneously compared with the fluid motions (see [48] for a suitable definition of T c ), a thin 
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reaction zone forms at the interface between burned and unburned material. It propagates 
into the surrounding fuel by one of two mechanisms allowed by the Rankine-Hugoniot 
jump conditions: a deflagration ("flame") or a detonation. 

If the overpressure created by the heat of the burning products is sufficiently high, a 
hydrodynamical shock wave forms that ignites the fuel by compressional heating. Such a 
self-sustaining combustion front that propagates by shock-heating is called a detonation. 
If, on the other hand, the initial overpressure is too weak, the temperature gradient at the 
fuel-ashes interface steepens until an equilibrium between heat diffusion (carried out pre- 
dominantly by electron-ion collisions) and energy generation is reached. The resulting 
combustion front consists of a diffusion zone that heats up the fuel to T c , followed by a 
thin reaction layer where the fuel is consumed and energy is generated. It is called a de- 
flagration or simply a flame and moves subsonically with respect to the unburned material 
[24]. Flames, unlike detonations, may therefore be strongly affected by turbulent velocity 
fluctuations of the fuel. Only if the unburned material is at rest, a unique laminar flame 
speed S\ can be found which depends on the detailed interaction of burning and diffusion 
within the flame region, e.g. [55]. Following [24], it can be estimated by assuming that in 
order for burning and diffusion to be in equilibrium, the respective time scales, Tb ~ e/w 
and Td ~ S 2 / k, where 6 is the flame thickness and k is the thermal diffusivity, must be 
similar: Tb ~ Td. Defining Si = 5/rb, one finds Si ~ (nw/e) 1 / 2 , where w should be 
evaluated at T w T c [48]. This is only a crude estimate due to the strong T-dependence of 
w. Numerical solutions of the full equations of hydrodynamics including nuclear energy 
generation and heat diffusion are needed to obtain more accurate values for S\ as a func- 
tion of p and fuel composition. Laminar thermonuclear carbon and oxygen flames at high 
to intermediate densities were investigated by [4, 14, 51], and, using a variety of different 
techniques and nuclear networks, by [48]. For the purpose of SN la explosion modeling, 
one needs to know the laminar flame speed Si w 10 7 . . . 10 4 cm s _1 for p w 10 9 . . . 10 7 g 
cm~ 3 , the flame thickness S = 10~ 4 ... 1 cm (defined here as the width of the thermal pre- 
heating layer ahead of the much thinner reaction front), and the density contrast between 
burned and unburned material p = Ap/p = 0.2 .. . 0.5 (all values quoted here assume a 
composition of Xc = Xq = 0.5 [48]). The thermal expansion parameter p reflects the 
partial lifting of electron degeneracy in the burning products, and is much lower than the 
typical value found in chemical, ideal gas systems [50]. 

Observed on scales much larger than 5, the internal reaction-diffusion structure can be 
neglected and the flame can be approximated as a density jump that propagates locally with 
the normal speed S\. This "thin flame" approximation allows a linear stability analysis of 
the front with respect to spatial perturbations. The result shows that thin flames are linearly 
unstable on all wavelengths. It was discovered first by Landau (1944) [23] and Darrieus 
(1938) [10] and is hence called the "Landau-Darrieus" (LD) instability. Subject to the LD 
instability, perturbations grow until a web of cellular structures forms and stabilizes the 
front at finite perturbation amplitudes [54]. The LD instability therefore does not, in gen- 
eral, lead to the production of turbulence. In the context of SN la models, the nonlinear 
LD instability was studied by [3], using a statistical approach based on the Frankel equa- 
tion, and by [33] employing 2D hydrodynamics and a one-step burning rate. Both groups 
concluded that the cellular stabilization mechanism precludes a strong acceleration of the 
burning front as a result of the LD instability. However, [3] mention the possible break- 
down of stabilization at low stellar densities (i.e., high p) which is also indicated by the 
lowest density run of [33]. 
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2.2. Hydrodynamic instabilities and turbulence 

The best studied and probably most important hydrodynamical effect for modeling SN 
la explosions is the Rayleigh-Taylor (RT) instability resulting from the buoyancy of hot, 
burned fluid with respect to the dense, unburned material [29, 30, 27, 18, 19, 34], and after 
more than five decades of experimental and numerical work, the basic phenomenology 
of nonlinear RT mixing is fairly well understood [11, 25, 45, 42, 53]: Subject to the RT 
instability, small surface perturbations grow until they form bubbles (or "mushrooms") 
that begin to float upward while spikes of dense fluid fall down. In the nonlinear regime, 
bubbles of various sizes interact and create a foamy RT mixing layer whose vertical extent 
/irt grows with time t according to a self-similar growth law, h^r = cigin/tyt 2 , where a 
is a dimensionless constant (a w 0.05) and g is the background gravitational acceleration. 

Secondary instabilities related to the velocity shear along the bubble surfaces [35] quickly 
lead to the production of turbulent velocity fluctuations that cascade from the size of the 
largest bubbles (rs 10 7 cm) down to the microscopic Kolmogorov scale, w 10~ 4 cm 
where they are dissipated [34, 19]. Since no computer is capable of resolving this range 
of scales, one has to resort to statistical or scaling approximations of those length scales 
that are not properly resolved. The most prominent scaling relation in turbulence research 
is Kolmogorov's law for the cascade of velocity fluctuations, stating that in the case of 
isotropy and statistical stationarity, the mean velocity v of turbulent eddies with size I 
scales as v ~ Z 1 / 3 [22]. Given the velocity of large eddies, e.g. from computer simula- 
tions, one can use this relation to extrapolate the eddy velocity distribution down to smaller 
scales under the assumption of isotropic, fully developed turbulence [34]. Knowledge of 
the eddy velocity as a function of length scale is important to classify the burning regime 
of the turbulent combustion front [37, 36, 20]. The ratio of the laminar flame speed and the 
turbulent velocity on the scale of the flame thickness, K = S\/v{8), plays an important 
role: if K > 1, the laminar flame structure is nearly unaffected by turbulent fluctuations. 
Turbulence does, however, wrinkle and deform the flame on scales I where S\ <C v(l), i.e. 
above the Gibson scale l g defined by Si — v(l g ) [40]. These wrinkles increase the flame 
surface area and therefore the total energy generation rate of the turbulent front [9]. In 
other words, the turbulent flame speed, S t , defined as the mean overall propagation veloc- 
ity of the turbulent flame front, becomes larger than the laminar speed S\. If the turbulence 
is sufficiently strong, v(L) ^> Si, the turbulent flame speed becomes independent of the 
laminar speed, and therefore of the microphysics of burning and diffusion, and scales only 
with the velocity of the largest turbulent eddy [9, 5]: 

S t ~ v(L) . (1) 

Because of the unperturbed laminar flame properties on very small scales, and the wrin- 
kling of the flame on large scales, the burning regime where K ^> 1 is called the corrugated 
flamelet regime [41, 5]. 

As the density of the white dwarf material declines and the laminar flamelets become 
slower and thicker, it is plausible that at some point turbulence significantly alters the ther- 
mal flame structure [20, 37]. This marks the end of the flamelet regime and the beginning 
of the distributed burning, or distributed reaction zone, regime, e.g. [41]. So far, modeling 
the distributed burning regime in exploding white dwarfs has not been attempted explicitely 
since neither nuclear burning and diffusion nor turbulent mixing can be properly described 
by simplified prescriptions. Phenomenologically, the laminar flame structure is believed 
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to be disrupted by turbulence and to form a distribution of reaction zones with various 
lengths and thicknesses. In order to find the critical density for the transition between both 
regimes, we need to formulate a specific criterion for flamelet breakdown. A criterion for 
the transition between both regimes is discussed in [37, 36] and [20]: 

'cutoff < 5 . (2) 

Inserting the results of [48] for S\ and S as functions of density, and using a typical turbu- 
lence velocity v(10 6 cm) <~ 10 7 cm s _1 , the transition from flamelet to distributed burning 
can be shown to occur at a density of pais ~ 10 7 g cm~ 3 [36]. 

3. A NUMERICAL SCHEME FOR TURBULENT COMBUSTION 

Our aim is now to convert some of the ideas presented in the previous Sections into 
a numerical scheme. The basic ingredients will be a finite-volume method to solve the 
fluid-dynamics equation, a front-tracking algorithm which allows us to propagate the ther- 
monuclear flame (assumed to be in the flamelet regime), and a model to determine the 
turbulent velocity fluctuations on unresolved sub-grid scales. 

3.1. The level set method 

The central aspect of our front tracking method is the association of the front geometry 
(a time-dependent set of points T) with an isoline of a so-called level set function G: 

T:={r\G(r) = 0} (3) 

Since G is not completely determined by this equation, we can additionally postulate that 
G be negative in the unburnt and positive in the burnt regions, and that G be a "smooth" 
function, which is convenient from a numerical point of view. This smoothness can be 
achieved, for example, by the additional constraint that | VG| = 1 in the whole computa- 
tional domain, with the exception of possible extrema and kinks of G. The ensemble of 
these conditions produces a G which is a signed distance function, i.e. the absolute value 
of G at any point equals the minimal front distance. The normal vector to the front is 
defined to point towards the unburnt material. 

The task is now to find an equation for the temporal evolution of G such that the zero 
level set of G behaves exactly as the flame. Such an expression can be obtained by the 
consideration that the total velocity of the front consists of two independent contributions: 
it is advected by the fluid motions at a speed v and it propagates normal to itself with a 
burning speed s. 

Since for deflagration waves a velocity jump usually occurs between the pre-front and 
post-front states, we must explicitly specify which state v and s refer to; traditionally, the 
values for the unburnt state are chosen. Therefore, one obtains for the total front motion 

Df — v u + s u n. (4) 

The total temporal derivative of G at a point P attached to the front must vanish, since G 
is, by definition, always at the front: 



(5) 
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FIG. 1. Illustration of the basic principles of the level set method according to [46]: The piecewise linear 
front cuts the mixed cells into burnt and unburnt parts, a is the unburnt volume fraction of a cell, /3 is the unburnt 
area fraction of a cell interface. The fluxes F u and F b are calculated from the reconstructed states. 



This leads to the desired differential equation describing the time evolution of G: 

3C 

- = -ZVVG. (6) 

This equation, however, cannot be applied on the whole computational domain, mainly 
because using this equation everywhere will in most cases destroy G's distance function 
property. Therefore additional measures must be taken in the regions away from the front 
to ensure a "well-behaved" |VG| [44]. 

The situation is further complicated by the fact that the quantities v u and s u which 
are needed to determine Df are not readily available in the cells cut by the front. In a 
finite volume context, these cells contain a mixture of pre- and post-front states instead. 
Nevertheless one can assume that the conserved quantities (mass, momentum and total 
energy) of the mixed state satisfy the following conditions: 

p = ap u + (1 - a)p b (7) 
pv = ap u v u + (1 - a)pbVb (8) 
pe = ap u e u + (1 - a)p b e h (9) 



Here a denotes the volume fraction of the cell occupied by the unburnt state. In order to 
reconstruct the states before and behind the flame, a nonlinear system consisting of the 
equations above, the Rankine-Hugoniot jump conditions and a burning rate law must be 
solved. 

Having obtained the reconstructed pre- and post-front states in the mixed cells, it is not 
only possible to determine Df, but also to separately calculate the fluxes of burnt and 
unburnt material over the cell interfaces. Consequently, the total flux over an interface can 
be expressed as a linear combination of burnt and unburnt fluxes weighted by the unburnt 
interface area fraction (3 (see Fig. 1): 

F = pF u + (l- p)F b . (10) 



3.2. Implementation 

For our calculations, the front tracking algorithm was implemented as an additional 
module for the hydrodynamics code PROMETHEUS [12]. Here we describe a simple 
implementation of most of the ideas discussed in the previous section which we will call 
"passive implementation". It assumes that the G-function is advected by the fluid motions 
and by burning and is only used to determine the source terms for the reactive Euler equa- 
tions. It must be noted that there exists no real discontinuity between fuel and ashes in 
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this case; the transition is smeared out over about three grid cells by the hydro-dynamical 
scheme, and the level set only indicates where the thin flame front should be. However, the 
numerical flame is still considerably thinner than in the reaction-diffusion approach. 

A complete implementation would contain in-cell-reconstruction and flux-splitting as 
proposed by [46] and outlined above. Therefore it should exactly describe the coupling 
between the flame and the hydrodynamic flow. This generalized version of the code has 
so-far only been applied to hydrogen combustion in air and, therefore, will be discussed 
elsewhere. 

3.2.1. G- Transport 

Since the front motion consists of two distinct contributions, it is appropriate to use an 
operator splitting approach for the time evolution of G. The advection term due to the fluid 
velocity vp reads in conservation form 

^S^d^r+l —vppGdf = (11) 
9t Jqv 

[28]. This equation is identical to the advection equation of a passive scalar, like the con- 
centration of an inert chemical species. Consequently, this contribution to the front prop- 
agation can be calculated by PROMETHEUS directly. The additional flame propagation 
due to burning is calculated at the end of each time step and a re-initialization of G is done 
in order to keep it a signed distance function (see [44]). 

3.2.2. Source terms 

After the update of the level set function in each time step, the change of chemical 
composition and total energy due to burning is calculated in the cells cut by the front. 
In order to obtain these values, the volume fraction a occupied by the unburnt material 
is determined in those cells by the following approach: from the value Gij and the two 
steepest gradients of G towards the front in x- and y-direction a first-order approximation 
G of the level set function is calculated; then the area fraction of cell ij where G < can 
be found easily. Based on these results, the new concentrations of fuel, ashes and energy 
are obtained: 

^Ashes = max ( 1 - a, X Ashes ) (12) 

Xpuel = 1 — -^Ashes (13) 
e[ ot = etot + ^(^Ashes ^ ^Ashes) (14) 

In principle this means that all fuel found behind the front is converted to ashes and the 
appropriate amount of energy is released. The maximum operator in eq. (12) ensures that 
no "reverse burning" (i.e. conversion from ashes to fuel) takes place in the cases where 
the average ash concentration is higher than the burnt volume fraction; such a situation 
can occur in a few rare cases because of unavoidable discretization errors of the numerical 
scheme. 



3.3. Turbulent nuclear burning 

The system of equations described so-far can be solved provided the normal velocity of 
the burning front is known everywhere and at all times. In our computations it is deter- 
mined according to a flame-brush model of Niemeyer and Hillebrandt [34], which we will 
briefly outline for convenience. 
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As was mentioned before, nuclear burning in degenerate dense matter is believed to 
propagate on microscopic scales as a conductive flame, wrinkled and stretched by local 
turbulence, but with essentially the laminar velocity. Due to the very high Reynolds num- 
bers macroscopic flows are highly turbulent and they interact with the flame, in principle 
down to the Kolmogorov scale. This means that all kinds of hydrodynamic instabilities 
feed energy into a turbulent cascade, including the buoyancy-driven Rayleigh-Taylor in- 
stability and the shear-driven Kelvin-Helmholtz instability. Consequently, the picture that 
emerges is more that of a "flame brush" spread over the entire turbulent regime rather than 
a wrinkled flame surface. For such a flame brush, the relevant minimum length scale is 
the so-called Gibson scale, defined as the lower bound for the curvature radius of flame 
wrinkles caused by turbulent stress. Thus, if the thermal diffusion scale is much smaller 
than the Gibson scale (which is the case for the physical conditions of interest here) small 
segments of the flame surface are unaffected by large scale turbulence and behave as un- 
perturbed laminar flames ("flamelets"). On the other hand side, since the Gibson scale 
is, at high densities, several orders of magnitude smaller than the integral scale set by the 
Rayleigh-Taylor eddies and many orders of magnitude larger than the thermal diffusion 
scale, both transport and burning times are determined by the eddy turnover times, and the 
effective velocity of the burning front is independent of the laminar burning velocity. 

A numerical realization of this general concept is presented in [34]. The basic assump- 
tion was that wherever one finds turbulence this turbulence is fully developed and ho- 
mogeneous, i.e. the turbulent velocity fluctuations on a length scale I are given by the 
Kolmogorov law v(l) = v(L)(l/ L) 1 ' 3 , where L is the integral scale, assumed to be equal 
to the Rayleigh-Taylor scale. Following the ideas outlined above, one can also assume that 
the thickness of the turbulent flame brush on the scale I is of the order of I itself. With these 
two assumptions and the definition of the Gibson scale one finds for Z g ib s ^5 I ^5 L ~ Art 

( I V /3 

v(l) ~ u t (l) ~ u^) — (15) 

V 'gibs / 

and d t {l) ~ I, where w(? g ibs) = "lam defines ? g it, s , ui am is the laminar burning speed and 
ut(l) is the turbulent flame velocity on the scale I. 

In a second step this model of turbulent combustion is coupled to our finite volume hydro 
scheme. Since in every finite volume scheme scales smaller than the grid size cannot be 
resolved, we express Z g ;b s in terms of the grid size A, the (unresolved) turbulent kinetic 
energy q, and the laminar burning velocity: 

/u 2 \ 3/2 

l gibs = A(^-\ • (16) 

Here q is determined from a sub-grid model [6, 34] and, finally, the effective turbulent 
velocity of the flame brush on scale A is given by 

ttt (A) = max(ttiam, v(A),vrt), (17) 
with v(A) = ^/2q and vrt oc \/gA, where g is the local gravitational acceleration. 

4. APPLICATION TO THE SUPERNOVA PROBLEM 

We have carried out numerical simulations in cylindrical rather than in spherical coordi- 
nates, mainly because it is much simpler to implement the level set on a Cartesian (r,z) grid. 
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FIG. 2. Snapshots of the temperature and the front geometry at 0.4s for model B5 of Reinecke et al. (1999a) 
(right figure) and turbulent velocity fluctuations on the grid scale (left panel). The position of of the front is 
indicated by the dotted curve. 
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FIG. 3. Snapshots of the temperature and the front geometry at 1.05s, taken for the low-resolution model 
(left figure) and high-resolution run, respectively. 

Moreover, the CFL condition is somewhat relaxed in comparison to spherical coordinates. 
The grid we used in most of our simulations maps the white dwarf onto 256 x 256 mesh 
points, equally spaced for the innermost 226x226 zones by A=1.5T0 6 cm, but increasing 
by 10% from zone to zone in the outer parts. The white dwarf, constructed in hydrostatic 
equilibrium for a realistic equation of state, has a central density of 2.9T0 9 g/cm 3 , a radius 
of 1.5T0 8 cm, and a mass of 2.8T0 33 g, identical to the one used by [34]. The initial mass 
fractions of C and O are chosen to be equal, and the total binding energy turns out to be 
5.4- 10 50 erg. At low densities (p < 10 7 g/cm 3 ), the burning velocity of the front is set equal 
to zero because the flame enters the distributed regime and our physical model is no longer 
valid. However, since in reality some matter may still burn the energy release obtained in 
the simulations is probably somewhat too low. 

Here we present only the results for two models, one in which nuclear burning was 
ignited off-center in a blob and a second one in which initially five blobs were burned as 
an initial condition, and refer to [43] for simulations with other initial conditions. 
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First, in Fig.2 a snapshot of the "five-blob" model B5 is shown at t = 0.4s. The leftt panel 
gives the position of the burning front, the temperatures (in gray-shading), and the expan- 
sion velocities. The rightpanel shows the distribution of turbulent velocity fluctuations. 
We find that, in accord with one's intuition, most of the turbulence is generated in a very 
thin layer near the front. Since in the limit of high turbulence intensity the nuclear flames 
propagate with the turbulent velocity it is obvious that this propagation velocity exceeds 
the laminar flame speed. However, for most of the initial conditions we have investigated 
this increase was not sufficient to unbind the star. In fact, model B5 of [43] was the only 
one of the set that did explode. 

Moreover, we show the results obtained with 3 times higher resolution for comparison. 
Fig. 3 gives a snapshot taken at 1.05s after ignition for the low-resolution run (left figure) 
and the high-resolution run, respectively. Although the increase in spatial resolution is 
only a factor of 3, the right panel shows clearly more structure. This is an important ef- 
fect because in the flamelet regime the rate of fuel consumption, in first order, increases 
proportional to surface area of the burning front. The net effect is that the low-resolution 
model stays bound at the end of the computations, whereas the better resolved model ex- 
plodes with an explosion energy of about 2xl0 50 erg. Fig. 3 also demonstrates that the 
level-set prescription allows to resolve the structure of the burning front down almost to 
the grid scale, thus avoiding artificial smearing of the front which is an inherent problem 
of front-capturing schemes. We want to stress that this gain of accuracy is not obtained at 
the expense of smaller CFL time steps because in our hybrid scheme the hydrodynamics is 
still done with cell-averaged quantities. 

5. CONCLUSIONS 

In this article, we have discussed the physics of thermonuclear combustion in degenerate 
dense matter of C+O white dwarfs. It was argued that not all relevant length-scales of 
this problem can be numerically resolved and, therefore, we have presented a numerical 
model to describe deflagration fronts with a reaction zone much thinner then the cells 
of the computational grid. The new approach was applied to the simulate thermonuclear 
supernova explosions. Our code works in 2 and 3 dimensions, but has only been used so 
far in 2D. 

An implicit assumption of our numerical model is that on the resolved scales the flows 
are turbulent allowing us to describe the physics on the unresolved scales by a sub-grid 
model, in the spirit of large-eddy simulations. The results presented here indicate that 
for supernova simulations this assumption is still not satisfied if (in 2D and cylindrical 
coordinates) we use a 256 x 256 grid. In fact, we could demonstrate that for a particular 
set of initial conditions increasing the numerical resolution by a factor of 3 only changes a 
model that remains bound after most of the nuclear fuel was burnt into a, although weak, 
explosion. The reason is simply that more structure on small length scales increases the 
rate of fuel consumption and we suspect that even our better resolved simulation does not 
yet give the final answer, pushing even 2D supernova simulation to the limits of what can 
be done on present supercomputers. 

The numerical scheme outlined in this paper can, of course, also be applied to chemical 
combustion in the flamelet regime and first successful attempts to model hydrogen com- 
bustion in air have already been made [46, 44]. Its advantage is always that it resolves the 
flame-front down to the grid-scale without making the time-steps intolerably small. 
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